Dynamical phase transition in correlated fermionic lattice systems 
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We use non-equilibrium dynamical mean-field theory to demonstrate the existence of a critical 
interaction in the real-time dynamics of the Hubbard model after an interaction quench. The critical 
point is characterized by fast thermalization and separates weak-coupling and strong-coupling 
regimes in which the relaxation is delayed due to prethermalization on intermediate timescales. 
This dynamical phase transition should be observable in experiments on trapped fermionic atoms. 
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The properties of correlated many-particle systems can 
change dramatically and abruptly as external parameters 
are varied. An important example is the Mott transition 
from an itinerant state to a correlation-induced insulator, 
which occurs in such diverse systems as transition-metal 
compounds [l| and ultracold quantum gases 0, S H]- 
An entirely new perspective on those systems is provided 
by their nonequlibrium dynamics after an external per- 
turbation, which is experimentally accessible not only in 
the case of well-controlled ultracold quantum gases, but 
also for electrons in solids by means of femtosecond spec- 
troscopy 0, 0, 0j B Hi ■ On short timescales the perturbed 
systems are essentially decoupled from the environment 
and follow the unitary time evolution according to the 
Schrodinger equation, which immediately raises a num- 
ber of questions: How does an isolated many-body sys- 
tem approach a new equilibrium after being quenched, 
i.e., after a sudden change in one of its parameters? Does 
it eventually thcrmalizc, or is detailed memory on the ini- 
tial state retained for all times? 

Recently these questions have been addressed in a 
number of experimental IC , and theoretical inves- 
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tigations [H, [13, 0, E 
a quench to a large interaction parameter U character- 
istic coUapsc-and-rcvival oscillations with period 2Tih/U 
appear, which are due to the integer eigenvalues of the 
interaction operator 0, [H, [li, [H, H 111, 0, [HI, [l!]. 
These oscillations eventually fade out, but in some cases 
the system is trapped in a nonthermal stationary state 
up to the largest accessible times 14 , l3| • At first glance 
this may seem surprising because thermalization is only 
known to be inhibited for integrable systems [l2, 13, 19| . 
whereas nonintcgrablc systems such as those studied in 
Refs. [3 15 are expected to thermalize 17 1. Nonther- 



mal quasistationary states are found, on the other hand, 
on an intermediate timescale (oc 1/C/^) after quenches to 
small interaction parameters [ist ; the prethermalization 
to these states is followed by an approach to thermal equi- 
librium on a much longer timescale (oc l/[/^). Another 
interesting observation is that the relaxation dynamics 
can depend sensitively on the parameters of the final 



Hamiltonian [2l|. Whether and how this phenomenon, 
which may be called a dynamical phase transition, re- 
lates to the existence of an equilibrium thermodynamic 
phase transition, remains to be clarified. 

Here we consider the relaxation of correlated lattice 
fermions described by a time-dependent Hubbard Hamil- 
tonian at half-filling. 
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using nonequilibrium dynamical mean-field theory 
(DMFT). We restrict ourselves to the paramagnetic 
phase and choose hoppings Vij corresponding to a semi- 
elliptic density of states p(e) = ^/W^~~e^ / {2iiV). The 
system is initially in the ground state of the nonintcract- 
ing Hamiltonian, i.e., J7(t<0) = 0. At t = the Coulomb 
repulsion is switched to a finite value, U (i>0) = U . En- 
ergy is measured in units of the quarter-bandwidth V 
and time in units of \/V, i.e., we set h = \ and in the 
figures also V = \. Our results confirm the prethermal- 
ization for quenches to U <^ V and indicate a second 
prethermalization regime for U V , for which wc pro- 
vide a general perturbativc argument. At C/^^" = 3.2F 
we observe a sharp crossover between the two regimes, 
suggestive of a dynamical phase transition in the above 
sense. 

Nonequilibrium DMFT. — In equilibrium DMFT, 
which becomes exact in the limit of infinite dimen- 
sions [23[, the self-energy is local and can be calcu- 
lated from a single-site impurity model subject to a self- 
consistency condition Nonequilibrium DMFT is a 
reformulation for Green functions on the Keldysh con- 
tour [25,, i26|] , which maps the lattice problem ^ onto a 
single-site problem described by the action 

5 = E jdtdt' ct{t)K{t,t')c,{t')+ j dth,,{t), (2) 

where h\oc{t) — U{t)[n-^{t)'~j\[ni{t)—}-\ is the local 
Hamiltonian at half-filling. For a system that is pre- 
pared in equilibrium with temperature T at times t 
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FIG. 1: Momentum distribution n{ek,t) for quenches from U 
— to U — 3 (left panel) and U = 5 (right panel). 

< the time-contour C is chosen to run from t = 
along the real axis to imax, back to 0, and finally to 
—i/T along the imaginary axis. The action S deter- 
mines the contour-ordered Green function Ga{t,t') = 
T:T[Tce-''^c^{t)c+{t')]/Z and the self-energy Y.{t,t'). For 
the semi-elliptic density of states the selfconsistency con- 
dition reduces to K{t,t') = V'^Ga{t,t') The single- 
site problem can be solved, for example, by means of 
real-time Monte Carlo techniques p?!. l28j. We use the 
weak-coup ling continuous-time Monte Carlo (CTQMC) 
algorithm [28|, which stochastically samples a diagram- 
matic expansion in powers of the interaction part h\oc 
and measures observables such as the local Green func- 
tion Gcr{t,t'). The weak-coupling method is highly suit- 
able for initial noninteracting states, because the imagi- 
nary branch of the contour does not enter the CTQMC 
calculation. This allows us to study initial states at zero 
temperature by transforming imaginary times to real fre- 
quencies. Furthermore, the parameters of the algorithm 
can be chosen such that only even orders contribute to 
Ga-{t,t') at half-filling, which reduces the sign problem. 
Computational details are deferred to a separate publi- 
cation. 

Many observables of the lattice model can be calcu- 
lated from the local Green function Ga{t, t') and the self- 
energy 5]CT(t,t')- The two-particle correlation function 
T„{t,t') = (Tcc+ (t)[n,s(i) - \]c^a{t')) is obtained from 
the contour convolution = G^ * So- and yields the 
double occupation d{t) — {ni'^{t)nii{t)) &i t = t' . (Local 
quantities do not depend on the site index i for the homo- 
geneous phase.) Solving the the lattice Dyson equation 
{idt + fi — e — Ecr) * Ggcr = 1 yields the momentum- 
resolved equal-time Green function Geait,t) and thus 
the momentum distribution n{€k,t) = {'^Xai^)'^kryi^)) ^ '"^ 
well as the kinetic energy per lattice site, Ekin(t)/L 
= 2 Jde p{e) n{e,t) e. The total energy E = Ekmit) + 
UL [d{t) — 1/4] must be conserved after the quench, and 
we use this as a check for the numerical solution. 

Results. — As depicted in Fig.[Tl the momentum distri- 
bution n(e, t) evolves from a step function in the initial 
state to a continuous function of e at large times. Re- 
markably, its discontinuity at e = 0, which marks the 
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FIG. 2: Fermi surface discontinuity An and double occupa- 
tion d{t) after quenches to U < 3 (left panels) and U > 3.5 
(right panels). Horizontal dotted lines in the upper left panel 
are at the quasistationary value Aristat ~ 2Z — 1 predicted 
in Ref. [3, with the T — Oguasiparticle weight Z taken from 
equilibrium DMFT data [30]. Horizontal arrows indicate cor- 
responding thermal values dth of the double occupation, ob- 
tained from equilibrium DMFT using QMC. Inset: thermal 
value dth and dmed, the average of the first maximum and the 
second minimum of d{t), which provides an estimate of the 
stationary value dstat; black dotted lines are the respective 
results from the strong-coupling expansion (see text). 

Fermi surface in the initial state, remains sharp while 
its height decays smoothly to zero. For a noninteract- 
ing initial state at half- filling, the discontinuity An{t) = 
n{0~ ,t) — ri(0+,t) can be expressed as 

An{t) = \Gl%Jt,0)\\ (3) 

where GZ^^t^O) = -*e(i)({c^,(0), c+^(i)}) is the re- 
tardcd component of the momentum-resolved Green 
function. This shows that the collapse of the disconti- 
nuity An is closely related to the decay of electron and 
hole excitations which are created at time < = at the 
Fermi surface. 

We now use An{t) and d{t) to characterize the re- 
laxation after the quench. As shown in detail below, 
these functions behave qualitatively different in the weak- 
coupling and strong-coupling regimes, separated by a 
very sharp crossover at U^^'^ = 3.2V. We test for 
thermalization by comparing to expectation values in a 
grand-canonical ensemble with the same total energy. 

Weak-coupling regime, U < U^'^"'. — For quenches to 
U < 3V (Fig. [l^,c) we find that the double occupation 
d{t) relaxes from its initial uncorrelated value d(0) = 
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(n|)o(nj,)o = 1/4 [291 almost to its thermal value dth, 
while the Fermi surface discontinuity An(t) remains fi- 
nite for times t < 5/V. This confirms the predictions by 
Moeckel and Kehrein [3l of a quasistationary state which 
is formed on timescales on the order of V/U'^ and has dgtat 
= dth + 0{U^/V^) and finite An^tat = 1-2Z. Here Z 
is the quasiparticle weight in equilibrium at zero temper- 
ature and interaction U. Their prediction was based on 
a perturbative flow equation analysis for U <^ V and it 
was argued that full thermalization occurs only on much 
longer timescales on the order oi V^/W^. At t = 2/V 
our numerical data agree very well with the predicted 
value of Aristat for U < IV. Note that even for quenches 
to larger f/, a prethermalization plateau remains visible 
in Fig. [2t at roughly this value, although the timescales 
V/U^ and V^/lf^ are no longer well separated. 

Strong- coupling regime, U > 11^'^"'. — For quenches to 
large U we observe collapse-and-revival oscillations with 
approximate frequency 27r/[/ both in d{t) and An{t) 
(Fig. (Jbjd). This phenomenon is well understood in the 
atomic limit {V = 0)jWhere the propagator e~*^* is ex- 
actly 27r/L/-pcriodic [2|. As expected, these oscillations 
arc damped for nonzero V: at least for small times, they 
fall off on timescales on the order of 1/V. Interestingly, 
the first few oscillations of d{t) are not centered around 
the thermal value dth (solid arrows in Fig. [JJd), which is 
instead located close to the first minimum of d{t) This 
suggests a prethermalization regime also for U ^ V, as 
discussed in Ref. uA, where oscillations in d{t) are damped 
to a nonthermal quasistationary value on the timescale 
1/V, while full thermalization can only happen on the 
longer timescale U/V'^. 

We now show that this prethermalization regime is 
a general feature of fermionic Hubbard-typc models at 
strong coupling and calculate the double occupation in 
the quasistationary state. We use the standard unitary 
transformation A = e~^Ae^ [sH for which the double 
occupation D = ^it'^U '-'^ dressed fermions ci^ 
is conserved, \H, D] = 0. After decomposing the hop- 
ping term K = Y^i^ia^ij-y l^yta^3<y^ ii^to parts 
that change the double occupation by p, i.e., = 
E.,.(^.ja/F)c+ c,.(l - n,,)n,-, = (A'_)+ and = K 
— — the leading order transformation is S — 

h.c. + 0{V^/U^). For 
3'^*De-''"*)o/L, we ob- 



iV/U)K+ + {V/UnK+,Ko] 
the double occupation, d{t) - 
tain 
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Re[e^'^ RitV)] +o(^, 



where R{tV) = (e'*^^° K+e-'*^'^°)o/ L and 4tat = 
d{0) + {2V/U)Rc{K+/L)o. The error 0{tV^/U^), which 
is due to omitted terms in the exponentials e^*^*, is ir- 
relevant in comparison to the leading terms iit ^ U /V'^. 
It remains to show that (i) the envelope function R{tV) 
of the oscillating term decays to zero for t ^ 1/V, and 
(ii) the quasistationary value dstat differs from the ther- 



mal value dth. (i) Inserting an eigenbasis KQ\m) = km\m) 
yields E(iF) =S„_„(|n)(TO|)oe'*^(*'-"'-'=") {n\K+\m). In 
this expression all oscillating terms dephase in the long- 
time average [itI . , so that only energy-diagonal terms 
contribute to the sum. But from [ii'oj-D] = it follows 
that Z? is a good quantum number of |n) so that (n|A'+ 
= 0, and thus R{tV) vanishes in the long-time limit (if 
it exists and if accidental degeneracies between sectors 
of different D are irrelevant). From Eq. ([4]) we there- 
fore conclude that d{t) equals dgtat for times 1/V <C t ^ 
U/V"^, up to corrections of order 0{V'^ /U"^). (ii) For the 
quasistationary value dstat we obtain 



rfstat =d(0) - Ad, 



(5a) 



Ad = - ^ jl^ict^Cj^in^s - njgf)a , (5b) 

ij(T 



which applies to arbitrary initial states. For noninteract- 
ing initial states the expectation value in this expression 
factorizes; in DMFT Eq. ([5b|) then evaluates to Ad = 
n{l — n/2){V /U){K / L)q, i.e., it is proportional to the 
kinetic energy in the initial state. For the thermal value 
dth we expand the free energy in V/T^:, because the ef- 
fective temperature is much larger than V after a 
quench toU ':$> V . At half-filling we obtain dth — d{Q) 4- 
{V/U){K/ L)q; for noninteracting initial states in DMFT 
we thus find that Ad = d(0) — dgtat = M(0) — dth]/2, i.e., 
at times 1/V ^ < <C U/V'^ the double occupation has 
relaxed only halfway towards dth- 

The strong-coupling predictions for the prethermaliza- 
tion plateau agree with our numerical results, for which 
the center of the first oscillation in d{t) approaches dgtat 
for large U (inset in Fig. [^b) . The scenario also applies 
to interaction quenches in the half-filled Falicov-Kimball 
model in DMFT ^ and the 1/r Hubbard chain 
although thermalization is inhibited in these models: in 
both models the long-time limit of d(t-^cx)) can be ob- 
tained exactly and indeed agrees with dstat for U V. 
For quenches to large U in the free 1/r chain (with band- 
width 2ttV) Eq. dSEI yields Ad = (y/C7)(l - 2n/3)7r. For 
the Falicov-Kimball model Ad is half as big as for the 
Hubbard model because only one spin species contributes 
to the kinetic energy in the initial state. 

Critical region, U w U^y" = 3.2V.— The characteris- 
tic collapse-and-revival oscillations of the strong-coupling 
regime disappear for quenches to U between 3.3T^ and 
3V, as is apparent from the Fermi surface discontinuity 
Ani at its first revival maximum (Fig. [3^). This change 
in the short-time dynamics reflects a change in the na- 
ture of single-particle excitations [Eq. ([3])] . It occurs also 
in equilibrium even at very high temperatures, because 
\Gl°J{t — i')p becomes oscillatory upon the transfer of 
spectral weight to the Hubbard subbands at ±C/. Ad- 
ditionally the prethermalization plateau at Arigtat disap- 
pears between 3V and 3.3T^, so that the system can relax 
rapidly after quenches to U values in this range: momen- 
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FIG. 3; Left panel: (7-dependence of the prethermalization 
plateau, Ajistat, and the first revival maximum, Am, sug- 
gesting a dynamical phase transition at (7^^" — 3.2V. Right 
panel: momentum distribution at time t — 1.6 after a quench 
toll = 3.3V; the momentum distribution has essentially ther- 
malized. The shaded area marks the critical region near f/^'*^" 
= 3.2V with very fast thermalization. 



turn distribution and double occupation reach their re- 
spective thermal values already before the first expected 
collapse-and-revival oscillation at time 2tt/U (Fig. [3Jd). 
This suggests that the prethermalization regimes at weak 
and strong Hubbard interaction are separated by a spe- 
cial point, which we estimate from our data to be lo- 
cated at C/^^" = 3.2V, where relaxation processes on all 
energy scales become relevant. This sharp crossover is 
quite remarkable in view of the fact that the effective 
temperature after the quench is much higher than the 
critical cndpoint of the Mott metal-insulator transition 
in equilibrium (T^ « 0.055^ whereas = 0.84^ 
for U ~ 3.3V), so that in equilibrium metallic and in- 
sulating phases could hardly be distinguished. A similar 



critical behavior was found for quenches in Heisenberg 
chains [2lj. Therefore one might speculate that such dy- 
namical phase transitions arc a generic property of the 
noncquilibrium dynamics of correlated systems; however, 
this issue requires further study. In particular it would be 
interesting to see whether there is also an abrupt change 
in the long-time behavior at the same ?7^^". 

Conclusion. — We determined the real-time dynamics 
of the Hubbard model after a quench from the nonin- 
teracting state using nonequilibrium DMFT, for which 
CTQMC 
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|28| | proves to be a very suitable impurity 
solver. This method allows to investigate the transient 
dynamics of correlated fermions in a variety of contexts, 
e.g., to describe experiments with cold atomic gases 
or pump-probe spectroscopy on correlated electrons [If 
For the Hubbard model we found that rapid thermaliza- 
tion occurs after quenches to [/ w C/^^"; in fact this is 
one of the few cases [H , 17 1 where thermalization can be 
numerically observed for a non-integrable model. On the 
other hand, for quenches to very small or very large in- 
teractions, the system becomes trapped in quasistation- 
ary states on intermediate timescales. The phenomena 
discussed in this paper are manifest in the momentum 
distribution and double occupation; both quantities are 
directly accessible in experiments with cold atomic gases. 
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